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The existence of dedicated neuronal modules such as those organized in the cerebral 
cortex, thalamus, basal ganglia, cerebellum, or spinal cord raises the question of how 
these functional modules are coordinated for appropriate motor behavior. Study of human 
locomotion offers an interesting field for addressing this central question. The coordination 
of the elevation of the 3 leg segments under a planar covariation rule (Borghese et al., 
1996) was recently modeled (Barliya et al., 2009) by phase-adjusted simple oscillators 
shedding new light on the understanding of the central pattern generator (CPG) processing 
relevant oscillation signals. We describe the use of a dynamic recurrent neural network 
(DRNN) mimicking the natural oscillatory behavior of human locomotion for reproducing 
the planar covariation rule in both legs at different walking speeds. Neural network 
learning was based on sinusoid signals integrating frequency and amplitude features 
of the first three harmonics of the sagittal elevation angles of the thigh, shank, and 
foot of each lower limb. We verified the biological plausibility of the neural networks. 
Best results were obtained with oscillations extracted from the first three harmonics in 
comparison to oscillations outside the harmonic frequency peaks. Physiological replication 
steadily increased with the number of neuronal units from 1 to 80, where similarity index 
reached 0.99. Analysis of synaptic weighting showed that the proportion of inhibitory 
connections consistently increased with the number of neuronal units in the DRNN. 
This emerging property in the artificial neural networks resonates with recent advances 
in neurophysiology of inhibitory neurons that are involved in central nervous system 
oscillatory activities. The main message of this study is that this type of DRNN may 
offer a useful model of physiological central pattern generator for gaining insights in basic 
research and developing clinical applications. 
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INTRODUCTION 

Neuronal modules profoundly influence many aspects of motor 
behavior. However, little is currently known about the con- 
trol mechanisms that allow the coordination of these modular 
entities. In this theoretical context, human locomotion is a chal- 
lenging movement because of the numerous neuroanatomical 
modules implicated in the different aspects of whole body move- 
ment, ranging from the cyclic propulsion of the limb to balance 
control. In spite of these intricate movement components and 
neuronal modules involved, it has been proposed that the over- 
all control can be realized by reducing the number of degrees 
of freedom of the system into global variables (Bernstein, 1967; 
Lacquaniti et al, 1999, 2002; Flash and Hochner, 2005; Latash, 
2008). The fact that the elevation angles of the three main lower 
limb segments are coordinated during gait within a covariation 



plane (Borghese et al., 1996), forming an elliptic loop corrobo- 
rates the idea that control of locomotion is also submitted to the 
general law of reducing variables (Barliya et al., 2009). This gen- 
eral law is also applicable for different walking speeds (Bianchi 
et al., 1998), for forward and backward directions (Grasso et al., 
1998), rectilinear or curvilinear walking paths (Courtine and 
Schieppati, 2004), walking with stilts (Dominici et al., 2009; Leurs 
et al., 2011), or with a transfemoral prosthesis walk (Leurs et al., 
2012), with different levels of body weight unloading (Ivanenko 
et al., 2002) and for running (Ivanenko et al., 2007). Notably, 
this inter-segmental coordination is not present in newly walk- 
ing toddlers (Cheron et al., 2001a,b; Ivanenko et al., 2005), but 
covariation planarity rapidly emerges over the first few days of 
independent walking experience, while the shape of the loop 
and plane orientation "mature" gradually over several years. 
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This evolution indicates that the attractor plane and the shape 
of the loop are neurophysiologically defined, rather than being 
imposed by biomechanical constraints (see Hicheur et al., 2006; 
and Ivanenko et al., 2008 for discussion). More recently, the 
developmental study of this complex behavior in a new born 
baby (Dominici et al., 2011) has permitted revisiting the con- 
cept of locomotor modules coding for the control of movement 
primitives. 

This modular approach raises the question of the dynamic 
coordination of modules in the context of oscillatory prop- 
erties of neuronal ensembles. Indeed, the dynamical structure 
of these modules must logically obey a common principle for 
movement generation: the production of oscillatory activity. 
Although this principle is well accepted in case of the rhyth- 
mic nature of locomotion (Georgopoulos and Grillner, 1989; 
Grillner, 2006), the recent study of Churchland et al. (2012) 
surprisingly demonstrates that non-periodic movements such 
as reaching are also generated by neuronal oscillation. This 
means that there is a strong possibility that the spinal mod- 
ules organized in a central pattern generator (CPG) could be 
dynamically controlled by cortical and/or supraspinal oscilla- 
tions. Interestingly, Barliya et al. (2009) recently modeled the 
time course of elevation angles of the three lower limb seg- 
ments in terms of simple oscillators coupled with appropriate 
time shifts reproducing the orientation of the plane and their 
elliptical shape. The oscillators were obtained by taking, after 
Fourier transform, the first three harmonics of the elevation 
angle kinematics. Each of these oscillators could be interpreted 
in term of oscillatory module acting in such a way that the 
orientation of the plane and the elliptic shape of the coordi- 
nation are conserved. It could thus be possible that oscillatory 
signals derived from these harmonics are used for activating CPG 
modules. 

The existence of CPG in the spinal cord in mammals has been 
proposed a century ago (Brown, 1911). In essence, it represents 
a group of neurons acting as a network to produce coordinated 
patterns of rhythmic activity. New evidence has shown the pres- 
ence of CPG in the spinal cord in humans (Calancie et al., 1994; 
Bussel et al, 1996). The characteristics of such CPG modules 
are their adaptability and robustness that lead to the production 
of different gait patterns adapted to their current environmen- 
tal context. For example, young infants (less than 1-year old) 
are already able to waUc over a split-belt treadmill with differ- 
ent types of coupling (Yang et al., 2005). Some of them were 
even able to walk in opposite directions. Mimicking physiology, 
the robotic and neuroscientific community developed artificial 
CPGs that are commonly used to animate robots from multi- 
legged insect-like to humanoids (for a review see Ijspeert, 2008). 
In their pioneering work in the cat, Patla et al. (1985) pro- 
posed an analytic model of limb locomotor pattern generator 
based on recorded muscle activity induced by electrical stimu- 
lation over the mesencephalic locomotor region (MLR) in the 
decerebrate cat. In this model, locomotor like patterns of six 
muscles resulted from six independent oscillators with dedicated 
parameters. These oscillators were reduced to simple sine and 
cosine functions fed by a tonic input. Since then, different meth- 
ods have been used from coupled non-linear oscillators (Ijspeert 



et al, 2007; Duvinage et al., 2011, 2012a), to highly detailed bio- 
physics of small groups of interconnected neurons (Hellgren et al., 
1992) and rhythm genesis of larger groups of neurons (Wadden 
et al., 1997) mainly in animal models. While human locomo- 
tion has often been reproduced computationally in the robotic 
field using equations that are numerically integrated (Righetti 
et al, 2005; Righetti and Ijspeert, 2006; Ceccato et al, 2009; 
Duvinage et al., 2011, 2012a), few methods involving neuron 
modeling for human gait generation have been studied so far. 
Among them Prentice and coauthors (1998, 2001) have success- 
fully transformed fundamental timing signals (sine and cosine 
inputs) into individual muscles activation bursts related to gait 
locomotion at different speeds using a feedforward neural net- 
work. Our group used the electromyographic (EMG) signals of 
the lower limb muscles as input for a dynamic recurrent neu- 
ral network (DRNN) producing as output the lower kinematics 
during locomotion (Cheron et al., 2003, 2012). However, the pos- 
sibility to produce the motion of the lower limb segments by 
means of oscillations derived from the three harmonics of the 
Fourier transform of walking kinematics has not yet been assessed 
by means of the same DRNN. We show here that after learn- 
ing based on different walking velocities, the DRNN is able to 
reproduce the lower limb kinematics of both legs. The DRNN can 
also generalize to the unlearned walking velocities. The analysis 
of the required network structure (e.g., number of units, distri- 
bution of time constant, and synaptic sign) provides a basis for 
the discussion about the constraints required for the elaboration 
of a CPG. 

METHODS 
EXPERIMENTAL SETUP 

Seven healthy subjects aged from 25 to 35 years (mean age: 
28 years) participated in this experiment. The protocol con- 
sisted of walking on a treadmill at 1 1 different velocities [from 
1 (0.28 m/s) to 6 km/h (1.67 m/s) stepped by 0.5 km/h (0.14 m/s)] 
leading to 1 1 trials per subject (total of 77 trials over all subjects). 
Whole body kinematics was recorded using Vicon system with six 
infrared Bonita cameras recording at 100 Hz during 40 s for each 
trial. The tracking consisted of recording 29 markers placed over 
the whole body. The position of the markers in an orthogonal ref- 
erence was computed using a custom biomechanical model. From 
the position of the markers of both legs, the kinematic (elevation) 
angles relative to the vertical plane of the laboratory have been 
calculated bilaterally for thighs, shanks, and feet. In this study 
we performed two experiments. The first one, "proof of con- 
cept," was done to ascertain the feasibility of a DRNN to learn 
elevation angles from pure sine waves. This part only includes 
supervised learning on a single pattern (Figure lA). The second 
experiment was performed in order to study the possibility to 
learn multiple patterns of walking (i.e., velocities) and to predict 
kinematics from unlearned patterns (Figure IB). After that, the 
DRNN structures were analyzed. 

DYNAMICAL RECURRENT NEURONAL NETWORK 

The DRNN is capable of modeling time-varying input-outputs 
and has varying weights as well as varying time constants for 
the artificial neurons (Pearlmutter, 1990). The adaptive time 
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FIGURE 1 I Realization of a CPG module based on experimental 
recordings of human wall<ing. (A) In the first step, the elevation angles of 
the thigh, shank, foot of the two legs are transformed in the time frequency 
domain using fast Fourier transform (FFT). Then the signal characteristics are 
computed for the first three harmonics and back-transformed into temporal 



space using sine waves formulation. These data are fed as input to the DRNN 
learning the transformation in kinematic signals (elevation angles). (B) A 
successful learning will permit the DRNN to predict kinematics based on sine 
wave inputs only. These output signals can produce biologically plausible 
walking patterns in a virtual reality avatar or an actual robotic exoskeleton. 



constants make it dynamic (Draye et al., 1996). The DRNN is 
governed by the following equations: 

T,^ = -y,+P{xO+I, (1) 
at 

where _F(fl) is the squashing function _F(fl) = (1 + e^'')^\y,- is the 
state or activation level of unit i, I„ is an external input (or bias), 
and X, is given by: 

X, = ^ WijYi (2) 
i 

which is the propagation equation of the network (x, is called 
the total or effective input of the neuron i, w,j is the synaptic 
weight between units i and;'). The time constants T,- will act like a 
relaxation process. In order to make the temporal behavior of the 
network explicit, an error function is defined as: 

ti 

E= j q(y(t),t)dt (3) 

fo 

where to and f i give the time interval during which the correction 
process occurs. The function q{y{t), t) is the cost function at time 
t which depends on the vector of the neuron activations y and on 
time. We then introduce new variables p,- (called adjoint variables) 
that will be determined by the following system of differential 
equations: 

^ = - - E (4) 

with boundary conditions p,(ti) = 0. After the introduction of 
these new variables, we can derive the learning equations: 

- = -jy,P'(x,)P/t (5) 
8£ 1 } dyi 

to 



The sinusoid signals derived from the Fast Fourier transform 
(FFT) kinematic data are sent as input to a DRNN (Figures 1, 2 
cf. Experiments 1, 2) that has to learn from these data to pro- 
duce the kinematics specified as elevation angles (Figure 1 A). 
Successful trainings were also used to produce kinematic patterns 
from unknown inputs (Figure IB) aiming to produce walking for 
multiple purposes, such as virtual avatars or robotic exoskeletons. 
The training is supervised, involving learning rule adaptations of 
synaptic weights and time constant of each unit (Draye et al., 
1995, 1996). A specific training procedure using Almeida algo- 
rithm was used to optimize learning performance (Cheron et al., 
2011). The DRNN involves a looping mechanism (fully con- 
nected structure) which enables this network to learn and store 
information (memory). This equips the network with the abil- 
ity to model complex situations with multiple influences. The 
DRNN was successfully used to map EMG signals into kine- 
matics during walking (Cheron et al., 2003), for the identifica- 
tion of the triphasic EMG pattern in subjects performing fast 
elbow flexion movements (Cheron et al, 2007) or to predict 
specific muscular activity in elite fencers compare to amateurs 
(Cheron et al, 2011). 

EXPERIMENT 1: PROOF OF CONCEPT 

DRNN computation has been used to prove that simple sine 
waves with specified temporal characteristics can be used as input 
to an artificial neural network to be transformed into elevation 
angles obtained from kinematic recordings during locomotion 
(Figures). 

Input: construction of sine waves 

As the learning phase of the DRNN is a time-consuming process, 
we had to select appropriate sample from the whole available data 
as input. Moreover, even if the kinematics during walking is rela- 
tively stable, it may vary too much to feed the DRNN during the 
learning phase. For these reasons, we extracted two consecutive 
gait cycles from the 40 s of experimental data recorded in each 
trial (Figure 2A, black curves). They were chosen so as to be rep- 
resentative in terms of frequency of the whole recording set. Then, 
to determine the kinematic characteristics of gait, we transformed 
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FIGURE 2 I (A) Elevation angles (in degree) of the three segments (thigh, 
shank, and foot) for each leg as a function of time for a subject walking 
at 3km/h (0.83 m/s). The whole pattern is presented for duration of 5 s. 
The black lines represent a pattern for two gait cycles used to determine 
the FFT characteristics. (B) The mean FFT transformation for six joints 
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for 40 s (in gray) and the two gait cycles (in black). Note that the two 
gait cycle patterns are selected so as to preserve the FFT characteristics 
in terms of amplitude and frequency (stars) for the overall pattern of 
40 s. These characteristics are used as parameters to generate sine 
waves as input of the DRNN. 
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FIGURE 3 I DRNN learning method for Experiment 1. For each subject, 
we trained 100 DRNNs to learn a pattern of kinematic [corresponding to 
a velocity of 3 km/h (3.83m/s)l from sine waves (Equations 7 and 8). 
Three sets of learning were defined as the input differed (SEA, SEB, 
SEC). The structures of DRNNs were modeled with 30 hidden units for 



each set of training. Each network iterated 10,000 times to change its 
synaptic weights and time constants. For each subject and structure we 
selected the network with the highest SI value. For each condition, the 
design of the network is then processed with 6 inputs, 30 hidden units, 
and 6 outputs. 



the data of the lower limb segments elevation angles into the time 
frequenq^ domain using the Matlab jft function to get the FFT 
power amplitude and their related frequency values of the first 
three harmonic peaks (Figure 2B). It has been shown previously 
that the first two Fourier harmonics accounted for approximately 
98% of the experimental variance of the thigh, shank, and foot 



angles (Bianchi et al, 1998). We decided to create sine waves 
based on the characteristics of the first three harmonics to ensure 
that the signal proposed as learning input to the DRNN contains 
enough information to match the desired output precisely. 

We extracted the values of the three frequencies (fl, /2, /3) 
corresponding to the maximal amplitudes (al, fl2, a3). It was 
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verified that /2 = 2/1 and /3 = 3/1. We then artificially pro- 
duced six sinusoidal signals using frequency values as parameters 
(Equations 7 and 8). 

yfi,i = sin(2jt x fi x t) (7) 
yfi,2 = sin(— (2Tt X fi X t)) (8) 
For! = 1, 2, 3 

These six sine waves thus correspond to the frequency charac- 
teristics of the kinematic patterns that will be further used as 
pattern to be learned. For the sake of clarity we called the set of 
original inputs set SEA (for set of Equations A). Additionally, we 
produced six new sine waves using the preceding computations 
(Equations 7 and 8) where /l' = fl + 0.25 Hz and six other sine 
waves where /I" =/l — 0.25 Hz, respectively called SEB (for set 
of Equations B) and SEC (for set of Equations C). Please note 
that in the latter two cases the original relations /2 = 2/1 and 
/3 = 3/1 were respected and hence /2' = 2fl + 0.50 Hz; /3' = 
3/1 + 0.75 Hz in the set SEB and /2" = 2/1 - 0.50 Hz; /3" = 
3/1 — 0.75 Hz in the set SEC. Three different input sets (SEA, 
SEB, SEC) were thus defined for learning. 

Pattern to learn: kinematic data 

The pattern to learned corresponds to the elevation angles of the 
right and left thigh, shank, foot segments in the two gait cycles 
of a 3km/h (0.83 m/s) walk, normalized between —1 and 1. The 
outputs were the same for SEA, SEB, and SEC whereas inputs 
differed. 

DRNN learning phase 

It was expected that the DRNN would learn how to transform the 
input to output thanks to its dynamical and recurrent structure 
of 30 hidden neurons. For each of the seven subjects, the network 
iterated 10,000 times. This process was repeated 100 times, lead- 
ing to 100 different networks. At the end of the learning phase, 
we selected for each subject the network for which the differ- 
ence between predicted and real kinematics was minimal. This 
computation was performed by calculating a similarity index (SI) 
(Bengoetxea et al, 2010) defined by the following equation: 

'^^Pl(t)p2(t)dt 

SI = ^ T (9) 

' r— r— 1 ^ 

where Tc is the period of the limit cycle, p\ and p2 being syn- 
chronized, i.e., the matching between both patterns is based on 
the maximum of each pattern. Note that if both functions are 
identical, SI = I. SI was calculated with the recorded pattern of 
elevation angles and the output of the DRNN. 

EXPERIMENT 2: MULTIPLE PAHERN LEARNING AND PREDICTION 

In this experiment sine waves were modulated in frequency 
and amplitude and transformed into kinematic data using 



multi-pattern training. The prediction of kinematic pattern from 
unknown sine wave pattern was also tested. 

Input 

As for the proof of concept methods, we extracted two 
gait cycles of each trial for multiple velocities (in km/h) 
(v = {1.5, 2.5, 3.5, 4.5, 5.5}). We transformed the kinematic 
data into the time -frequency domain to get their frequency 
(fl, fl, /3) and amplitude (al, al, a?!) (Figure 2B) parameters 
using the following set of Equations (10 and 11). 

yy,fiv, fliv.i = aiy x (sm(2jt x fi^ x t)) (10) 
yv,fiy, flZv,2 = aiy X (sOT(-(2n x /^ x f))) (H) 
For! = 1, 2, 3 

Patterns to learn: kinematic data 

The patterns to be learned consisted of elevation angles of the 
right and left thigh, shank, and foot segments corresponding to 
the two gait cycles, normalized between -1 and 1. These patterns 
were obtained from recordings at multiple velocities (in km/h) 
(v= {1.5, 2.5, 3.5, 4.5, 5.5}). 

DRNN training phase 

We used the possibility to train the DRNN in a multi-pattern pur- 
pose (Bengoetxea et al, 2005). For each subject, the DRNN was 
trained to match the inputs/outputs patterns corresponding to 
five different velocities within a single DRNN structure. Hundred 
networks were made per subject, each of them iterating 50,000 
times. This operation was assigned for 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 
30, 40, 50, 60, 70, and 80 hidden units. When this was completed, 
we selected the best network for each subject and each number of 
hidden units using the maximal SI values. 

Output prediction 

We built sine waves from intermediary pattern velocities (km/h) 
(v = {2, 3, 4, 5}) as explained above. We fed the best DRNN 
structures previously obtained after the training phases with these 
unlearned inputs sine waves and analyzed DRNN behavior by 
calculating the predicted output with experimental data using SI 
values. 

NEURONAL PROPERTIES AND CONNECTIVITY AFTER LEARNING 

Introduction of timing allows modeling of more complex fre- 
quency behavior, improves the non-linearity effect of the sig- 
moid function and the memory effect of time delays (Draye 
et al, 1995). The distribution of the time constant and the 
synaptic weights between units (Draye et al., 1996) after learn- 
ing was analyzed after multiple pattern learning and predic- 
tion. In particular, we recorded the number of positive and 
negative connections inside DRNN structures of the best net- 
works. Additionally we studied the distribution of neuronal 
time-constants. 

STATISTICAL ANALYSIS 

Statistical analysis was performed using Statistica software 
(StatSoft, www.statsoft.com). For each test performed and 
described in the result section, we firstly verify whether the 
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data were distributed normally using Kolmogorov-Smirnov test. 
When the data were normally distributed we use ANOVA with 
repeated measures and post-hoc Fisher analyses. When it was not 
possible to use ANOVA we used non-parametric tests such as 
Friedman ANOVA or sign tests. 

COMPUTATIONS 

AH DRNN computations were performed in the Hydra com- 
puting center hosted in ULB (https://cc.ulb.ac.be/hydra/). We 
allocated 1 node and 10 Gb of memory per computation (i.e., per 
subject per condition in the Experiment 1, per subject per struc- 
ture in the Experiment 2). The computations were run in parallel 
to optimize the learning duration. Afterwards we simulated 5% 
of the overall experiment in the same conditions to estimate the 
learning time. The overall duration of the process was obtained 
by linear interpolation (Figure 4). 

RESULTS 

EXPERIMENT 1: PROOF OF CONCEPT 

A statistical test was designed to compare SI values from different 
input types (SEA, SEB, or SEC) (Figures 3,5). The Kolmogorov- 
Smirnov did not reject the hypothesis that SI values were 
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FIGURE 4 I Illustration of the DRNN structure and computational 
time. (A) Computational learning duration for one subject for 
experiment 2 (in hours). The gray circles represent the estimation of 
time required per number of hidden units from a recalculation of 5% 
of the procedure. Please note that the overall process for the seven 
subjects has taken approximately 200 days to be computed. (B), (C), 



normally distributed (D = 1.4414, p > 0.20) when analyzing 
together values of the different inputs (SEA, SEB, or SEC). We 
used an ANOVA with repeated measures where dependant vari- 
able was SI and the independent variable chosen was the type 
of input. The analysis showed an effect of the input frequency 
in the DRNN prediction (SI value) [H F(2, 12) = 38.110, p = 
0.00001]. Post-hoc analysis confirmed that SI values of original 
group of unchanged frequency input (SEA) were higher than the 
2 modified groups where frequency inputs have changed (SEB 
and SEC). 

EXPERIMENT 2: PREDICTION OF INTERMEDIARY VELOCITIES 

As we decided to use frequency (f l,/2,/3) and amplitude {al, al, 
fl3) characteristics to modulate inputs for multiple learning pro- 
cedures (Equations 10 and 11), we have verified that there was a 
statistical significance of these parameters for different velocities. 
The Kolmogorov-Smirnov test for /I (D = 0.14370, p < 0.1), 
/2 (D = 0.14370, p < 0.1), and /3 (D = 0.14370, p < 0.1) was 
not clear enough to reject the fact their population may fol- 
low a normal law. When looking at the distribution, they tend 
to be normal and it is possible that the significance of the test 
is due to the weak number of values. According to similar test. 



Hidden units connectivity 

and (D) represent the connectivity of the hidden layer with 5, 10, and 
20 neurons, respectively. As the number of neurons increases, the 
number of connections increases by a factorial multiplication as well as 
the time required to adjust connection weight and time constant. 
Please note that only the units of the hidden layers are represented 
without input or output neurons. 
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FIGURE 5 I Planar covariation between normalized tliigli, shanl<, and 
foot for all participants. (A) Recording data — the patterns correspond to the 
real kinematics of 2 gait cycles for all participants. Best kinematic patterns 



the values for parameters al {D = 0.08476, p > 0.2), a2 (D = 
0.08758, p > 0.2), and a3 (D = 0.09035, p > 0.2) were normally 
distributed. We then used two ANOVA with repeated measures 
where the dependent variables were, respectively, the amplitude 
and frequency values and the within-subject factors were velocity 
of walking and the specific harmonic (1, 2, or 3). 

ANOVA shows significant changes in amplitude [-F(io, 60) = 
31.351, p < 0.0001] and frequency [_F(io, eo) = 69.276, p < 
0.0001] with an increase in velocity. Post-hoc analyses revealed 
an increase in /I, /2, al, al, a3 and a decrease of/3 with an 
increase in the velocity. These significant differences justified their 
use for building specific sine waves for different walking velocities 
(Figure 6). 

Concerning the performance of the DRNN outputs, we 
wanted to verify if the SI value applied for the best networks was 
different for learning and prediction. Additionally we wanted to 
statistically check if the number of hidden units of the networks 
increases the SI value (Figure 7). The Kolmogorov-Smirnov did 



B 




-1 -1 



predicted by the DRNN for each subject where SEA (B), SEB (C), and SEC 
(D) were fed as input. SI is the average of similarity index value {SI) in each 
condition for all subjects when compared with the recording data. 



not verify that the populations of SI values among the learn- 
ing (D = 0.18467, p < 0.01) or the prediction (D = 0.20684, 
p < 0.01) were normal. Thus, to compare the SI values between 
learning and prediction process we chose to use a sign test as 
the structure of the network (weights and time constant) were 
the same. This test reveals no significant differences in SI values 
between the two populations except when the network contains 
1, 2, 4, 5, 6, and 7 neurons. Moreover we use a Friedman test to 
analyze SI values (dependent variable) according to the number 
of hidden units (independent variables) of the network. There is 
an effect of the number of hidden units to the SI both for learn- 
ing [x^g) = 111.5630, p < 0.0001] and for prediction [Xjig) = 
109.3721, p < 0.0001]. 

We also analyze the output of the DRNN for specific veloc- 
ities where the hidden layer was the biggest with 80 units. A 
Kolmogorov-Smirnov run together with learning and prediction 
values did not verify that the populations of SI was normal (D = 
0.18342,p < 0.05). 
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FIGURE 6 1 Evolution of the frequency (A) and amplitude (B) of 
the first three harmonics when transforming elevation angles 
into the time frequency domain (FFT) at different walking 
velocities [from 1 to 6 km/h (0.28 to 1.67 m/s)]. The FFT peaks 



of each angle are averaged for each velocity for both frequency 

and amplitude. Squares, circles, and triangles represent the mean 

value for 7 subjects and whiskers represent the 95% confidence 
interval. 



We use a Friedman test to analyze SI values (dependent vari- 
able) according to the velocity (independent variables) of the 
network. There is an effect of the velocity to the SI Ixls) ~ 30-22, 
p = 0.00019] (Figures). A post-hoc analysis at 0.05 significance 
level reveals that SI values of 4.5 km/h were different from 1 .5 and 
2 km/h. SI values of 4 km/h were also different from 2 km/h. An 
example of prediction of intermediate velocities in one subject is 
illustrated in (Figure 9). 

DRNN STRUCTURES FROM EXPERIMENT 2 
Weight distribution analysis 

The percentage of positive and negative weights was calculated 
for each best network per subject per condition. We wanted to 



verify if the distribution of the positive and negative weights 
were different (Figure lOA). The Kolmogorov-Smirnov test did 
not verify that the distribution of positive weights were normal 
(D = 0.15022, p < 0.01) nor the distribution of negative weights 
(D = 0.15022, p < 0.01). We then used a non-parametric sign- 
test to compare the two-distribution as the two samples were 
dependent. Regardless of the number of hidden units, the test 
shows a difference between the distributions of the two popu- 
lations (Z = 9.482, p < 0.0001). When the number of units in 
the test were included, it appears that the population of nega- 
tive weights is higher than population of positive weights when 
the network possesses more than 3 hidden units (for 4 hid- 
den units, Z = 2.268, p = 0.023342). When the structure of the 
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FIGURE 7 I SI values of all multi-patterns learning and prediction for all participants with different number of hidden neurons. Whiskers correspond to 
95% confidence interval. 
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FIGURE 8 I SI values for learning and prediction with a 80 hidden units network according to wall<ing velocity. Whiskers correspond to 95% 
confidence interval. 



DRNN reaches 80 hidden neurons, 70.6 ± 0.84% of the synapses 
are negative. 

Time constant distribution analysis 

The distribution of the time constant was represented by 
the median of neuronal time constants of the best networks 



per subject and number of hidden units (Figure lOB). The 
Kolmogorov-Smirnov displayed a normal distribution of the 
time constant median (D = 0.09110, p > 0.20). ANOVA with 
repeated measures was designed with the number of hid- 
den units as independent variable and time constant median 
as dependent variable. ANOVA analysis shows an interesting 
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FIGURE 9 I Prediction of intermediate velocities in one subject. 

(A) Planar covariation of the normalized thigh, shank, and foot for one 
participant walking at 2km/h. (B) Planar covariation of the normalized 
thigh, shank, and foot for the same participant walking at 5km/h. (C) 
Kinematics of walking of the normalized thigh, shank, and foot for the 



same participant walking at 2km/h. (D) Kinematics of walking of the 
normalized thigh, shank, and foot for the same participant walking at 
5km/h. In each subplot the red dotted line corresponds to the 
experimental data whereas black thick line corresponds to the output of 
the DRNN. 



effect of the number of hidden neurons [F(i6_ 95) = 3.6245, 
p < 0.0001]. Overall, post-hoc analysis revealed that the dis- 
tribution of the time constants were different when the net- 
work was small (less than 5 hidden units). It also reveals 
that the networks with 80 hidden neurons were differ- 
ently distributed than medium-sized network (8-10 hidden 
neurons). 

DISCUSSION 
MAIN FINDING 

We show here that a fuUy connected recurrent neural network 
is able to reproduce human walking pattern based on oscilla- 
tory properties of kinematics. Although, this network is a black 
box model without prewired structure mimicking a physiologi- 
cal CPG, its actual performance allows direct comparison with 
CPG dedicated structure and related algorithm (Duvinage et al., 
2012a). Moreover, by the inherent input-output mapping, the 



DRNN models not only the CPG but also neural feedback 
pathways and the musculoskeletal system. For simplicity, we con- 
sider this neural network as "CPG-like structure" here. We proved 
that the DRNN is capable of generating the kinematics as ele- 
vation angles pattern of walking for both limbs (six degrees of 
freedom) from simple oscillations corresponding to the three 
main harmonics of the walking kinematics. Moreover, by mod- 
ulating those frequencies and tuning them in amplitude as input, 
the DRNN was able to learn and reliably predict walking kine- 
matics at different velocities (Figure 9). After this appropriate 
learning the DRNN can thus be considered as a CPG-like struc- 
ture that would continuously receive oscillation inputs to produce 
the relevant elevation patterns of the six leg segments. Another 
interesting result is observed when looking at the structure of 
best CPGs obtained after learning. Hence, it appears that all of 
them contain a major part of negative connection weight between 
units. 
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FIGURE 10 1 Distribution of weights (A) and time constants (B) 
among the DRRN structures. The horizontal axis represents the 
number of hidden neurons and vertical the percentage of 



population. The bars represent the mean distribution among the 
best DRNN for each subject. Whiskers represent the 95% 
confidence intervals. 



LIMITATIONS OF THE PRESENT APPROACH 

Obviously, there are infinite different ways to train the DRNN 
and this is a strength of this approach. However, it implies a 
corollary limitation as not all possibilities could be tested in 
the present study. For example, in order to better document its 
generalization ability, the model could be trained for a defined 
low range of velocities, e.g., from 1 to 4 km/h, and then tested 
with unlearned oscillation input corresponding to higher veloc- 
ities. A reverse procedure could also be made, i.e., from faster 
training to slower predictions. Furthermore, inter- subjects gen- 
eralization has not been studied in the present investigation. 
The actual usefulness of performing this would largely depend 
on the basic or application purposes. Another limitation of the 
present work is the lack of feedback testing, which necessitates 
a priori identification of a reliable signal and a new operational 
strategy for learning. Future work will address these aspects 
specifically. 

NEUROPHYSIOLOGICAL SIMILARITY BETWEEN MODELED CPG AND 
CPG IN HUMANS AND OTHER MAMMALS 

The understanding of CPG mechanisms remains central in loco- 
motion study (Grillner, 2006; Kiehn, 2006; Rossignol and Frigon, 
2011). The CPG is a spinal network of neurons capable of gen- 
erating a rhythmic pattern of alternate activities between flexor 
and extensor motoneurons on the same side with reciprocal acti- 
vation of homologous motoneurons in the contralateral limb. 
This intrinsic spinal circuitry has been well described in many 
invertebrate and vertebrate animals, and is highly conserved even 
in humans, where greater cortical control of spinal modules is 
required working in conjunction with sensory feedback (Calancie 
et al., 1994; Bussel et al., 1996; Duysens and Van de Crommert, 
1998; Drew et al, 2002; Rossignol et al., 2009). The unique char- 
acteristics of human walking probably reflect a complex neural 
mechanism responsible for pattern production. It is therefore 
difficult to directly extend experimental findings obtained in 



quadruped animals to human walking (Barbeau et al, 1998; 
Capaday, 2002). 

The fact that some patients with incomplete spinal injury can 
move their legs in a rhythmic fashion (Dietz et al., 1995) and 
that the primary sensorimotor cortex provides oscillatory com- 
mands toward the spine during walking (La Fougere et al., 2010) 
motivates new experiments in which different types of oscillatory 
signals could be used as input to the CPG-like DRNN. In this con- 
text, recent studies have showed EEC oscillations in relation to 
the gait cycle phase including event-related spectral perturbation 
in the alpha-beta and gamma bands (Gwin et al., 2011; Haefeli 
et al, 2011; Cheron et al, 2012; Wagner et al, 2012). 

These results are consistent with a top down control of 
locomotion (Capaday, 2002) and demonstrate the feasibility of 
extracting EEG signals from the sensorimotor cortex controlling 
the contralateral foot placement during walking. Although the 
distinction between the brain signals directly linked to the motor 
commands and those related to the treatment of multiple sen- 
sory signals is a hard task. In this context, Petersen et al. (2012) 
have found evidence of synchrony in the frequency domain 
between the primary motor cortex and the tibialis anterior muscle 
prior to heel strike during the swing phase of walking signifying 
that rhythmical cortical activity is transmitted via the corti- 
cospinal tract to the active muscles. Additionally Wagner et al. 
(2012), showed a significant difference in the alpha (8-12 Hz) 
and beta (18-21 Hz) rhythm recorded over the central midline 
area between passive and active walking with exoskeleton. The 
role played by specific oscillations related to the initiation and 
control of human locomotion coming from supraspinal struc- 
ture was recently demonstrated by local field potential recordings 
performed in the pedunculopontine nucleus in parkinsonian 
patients during rest and unconstrained walking (Thevathasan 
et al., 2012). Alpha oscillation recorded in the caudal part of 
this nucleus is correlated with gait speed and permits to suppress 
"task irrelevant" distraction for improving gait performance. 
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Moreover, these authors showed that gait freezing of parkinsonian 
patient was associated with the attenuation of these alpha waves. 

Consistent with this aspect of gait physiology, in our model, 
input sine waves are sufficient to predict successful output with 
the DRNN and offer the possibility to mimic such type of 
supraspinal oscillatory input. It could be advanced that nonlinear 
mapping of sinusoidal oscillations to kinematic pattern should 
be realized by other mathematical functions, such as by a Taylor 
series, but such multi-dimensional approximation seems to be 
highly difficult to obtain and does not permit testing different 
network configurations mimicking biological organization such 
as the CPG. In the present study, we focused on input oscilla- 
tion derived from the first three harmonics of kinematic signals, 
which were slower than the alpha frequency range. However, it 
could be possible to extract slower oscillation from alpha or beta 
derived signals (envelope) in order to activate the DRNN in our 
future work. 

STRUCTURAL SIMILARITIES OF MODELED CPG AND 
NEUROPHYSIOLOGICAL CPGs IN ANIMALS 

It has been suggested that the biological CPG is considered to 
serve two basic functions: rhythm generation (RG) and pattern 
formation (PF). 

Initially proposed by Perret and Cabelguen (1980), the idea 
that the biological CPG is composed of a rhythm and a pattern- 
amplitude generator is now widely accepted (Kriellaars et al., 
1994; Guertin et al, 1995; Perreauk et al, 1995; Grillner, 2006; 
Kiehn, 2006; Talpalar et al., 2011) and paved the way to more 
complex models of multi-level CPG (McCrea and Rybak, 2008, 
see below). 

It is well-recognized that rhythm generating networks can be 
realized by means of (1) pacemaker neurons with intrinsic mem- 
brane properties such as those described in the stomatogastric 
ganglion of crustaceans or in the mammalian thalamus (Steriade 
and Llinas, 1988) or (2) most simplistic neurons without intrin- 
sic pacemaker properties but interacting with inhibitory synapses 
for producing oscillation as emergent properties of this neuronal 
population (Geisler et al., 2005). Both neuronal systems thus 
present the fundamental ability to oscillate. Firstly described in 
the tadpole and lamprey CPGs, glutamatergic excitatory neurons 
distributed along the cord (Grillner, 2003) assume the func- 
tion of rhythmic generator by driving motor neurons and other 
ipsilateral and commissural inhibitory neurons coordinating the 
different CPG modules. By blocking the inhibitory networking in 
the lamprey and also in rodent and cat, many authors (see Kiehn, 
2006 for a review) have demonstrated that the glutamatergic burst 
neurons are the generators of the CPG rhythm. 

In addition to intrinsic RG properties, the walking CPGs need 
to integrate the ipsilateral coordination of flexors and extensors 
across the same or different joints in a limb and perform inter- 
limbs coordination. It has been proposed (Zhong et al., 2012) that 
a subpopulation of neuronal CPG that drives extensor activity is 
tonically active and is regulated via inhibitory interactions with 
another CPG rhythmic structure responsible for flexors activity 
in the same hemicord. This assumption may explain why, during 
experimental recordings on the neonatal mouse isolated spinal 
cord, spontaneous deletions of extensor activity do not perturb 



rhythmic flexor activity. Thus, the inhibitory interneurons play 
a major role in the temporal sculpting and coordination of the 
CPG units. The interneurons and the Renshaw cells are involved 
in this function and in the regulation of walking speed. In addi- 
tion, the left-right coordination is assumed by a complex network 
of excitatory and inhibitory commissural interneurons acting on 
both motor neurons and inhibitory interneurons of the con- 
tralateral side (Kiehn, 2006). Interestingly, we have shown that 
a great percentage of artificial neurons became inhibitory neu- 
rons (negative synaptic weight) when the number of neurons 
progressively increases in the DRNN structure. In this context, 
it was recently demonstrated in awake mice that the spiking 
activity of inhibitory neurons of the barrel cortex is organized 
in order to balance excitation and prevent explosive activity in 
the recurrently connected cortical microcircuit (Gentet et al, 
2010). This physiological mechanism can also be proposed in 
the present case of the emergent structure of artificial DRNN 
circuit. Another, not exclusive explanation can reside in the preva- 
lence of inhibitory recurrent connections for producing network 
oscillation (Geisler et al., 2005; Wildie and Shanahan, 2011). 
In their review, Nishimaru and Kakizaki (2009) have proposed 
that inhibitory interneurons play a major role in the CPG of 
rodent spinal cord. The interneurons are likely to control the 
bursting of motor neurons during locomotion and it appears 
that the synaptic transmission mediated by glycine and GABA 
shifts from excitatory to inhibitory during the prenatal period. 
It was recently demonstrated that in the absence of glutamater- 
gic synaptic transmission, the flexor-extensor alternation appears 
to be generated by the inhibitory interneurons, mediating recip- 
rocal inhibition from muscle proprioceptors to antagonist motor 
neurons (Talpalar et al, 2011). 

The present artificial model does not pretend to mimic the 
complexity of the CPG structure. Instead, it presents a highly 
simplified recurrent organization from which CPG-like dynamic 
function emerges, following appropriate learning. Sinusoidal 
inputs serve as temporal referent to produce rhythmic angles pat- 
terns. This model could correspond to the RG structure described 
previously as a higher order structure that determines rhythmic 
output of the system (McCrea and Rybak, 2008; Zhong et al, 
2012) since sine waves are transformed into kinematics. Another, 
lower order structure responsible for the phasing and intensity 
coordination (McCrea and Rybak, 2008; Zhong et al., 2012) could 
be assumed by another model of DRNN transforming theoretical 
kinematics into practical muscle orders. We have already stud- 
ied such relation where EMG signals from walking where used 
to predict kinematics (Cheron et al., 2003). To conclude to this 
point, we propose that the two model driving two specific DRNNs 
(one for producing elevation angles from sine waves and one pro- 
ducing muscular patterns from elevation angles during walking) 
could act as a complementary top-down pathway to produce ade- 
quate coordinated patterns as it has been proposed to model the 
locomotion in spinal mouse (Zhong et al, 2012). 

The present results can also be discussed in the light of the 
electrical stimulations performed in the MLR) inducing locomo- 
tor behavior in decerebrate cat (Shik et al., 1969) or in lamprey 
(McClellan and Grillner, 1984). In mice, prolonged rhythmic 
stimulation on the midline of the caudal hindbrain or the ventral 
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spinal cord (C1-C4) induces a stable locomotor activity (Talpalar 
et al., 2011). Typically, low-frequency stimulation leads to slow- 
frequency movements and inversely fast-frequency stimulation 
leads to fast- frequency movements. Our model is in accordance 
with this physiological behavior, when amplitude of the artificial 
sine wave inputs increases, the amplitude of stepping increases 
as well, leading to a change in walking velocities (Figure 6B). In 
terms of neurological development, there is some evidence for the 
existence of CPG very early in CNS maturation (Yang et al, 1998; 
Dominici et al., 2011). Neonatal, so-called "infant" stepping has 
been ascribed to similar EMG patterns activity in different direc- 
tions inducing stereotyped yet non-functional walking patterns 
(Lamb and Yang, 2000). This leads the authors to conclude that 
the same CPG controls different stepping in human infants in 
contrast with some studies in adults (Thorstensson, 1986; Grasso 
et al, 1998). Interestingly, we found that DRNN with only four 
hidden artificial neurons can generate walking pattern, whereas 
at least 50 hidden neurons are the prerequisite to generate accu- 
rate movements (Figure?). Obviously, recruitment and training 
of such high numbers of neurons requires long computational 
time. For example, with a 4 hidden units DRNN structure, the 
learning process lasts about 5 min based on our computer per- 
formance while 160 min are necessary for a DRNN containing 
80 units (Figure 4). 

PERSPECTIVES 

Such a tool can be used to produce gait kinematics in numer- 
ous and various applications. For rehabilitation it can be used 
to train people for recovering a walking pattern corresponding 
to their physical characteristics by training with an appropri- 
ate feedback. Specifically dedicated DRNN based on the proper 
dynamics of participants could be used for medical applications 
such as in prosthesis and exoskeleton control (Cheron et al., 
2012). It can also be integrated in BCI applications in which 
higher order commands can be used, e.g., from steady state visual 



or somatosensory evoked potentials (Cheron et al, 2012) or P300 
(Castermans et al., 2011; Duvinage et al., 2012b). This neuronal 
avenue might lead to the decoding of higher neuronal commands 
that govern CPG mechanisms. Since these CPG can be trained 
using specific sinusoidal frequency signals, it might be possible 
to extract this type of signals from specific EEC rhythms since 
the brain itself is an effective machine for producing oscilla- 
tions (Buzsaki and Draguhn, 2004). One of the strengths of this 
approach is that it is not necessary to determine in advance the 
topology and the timing sequences between the artificial neu- 
rons. This contrasts with other CPGs, such as a recently developed 
ones (Duvinage et al., 201 1) based on coupled oscillators (Righetti 
and Ijspeert, 2006), where adjustment of intrinsic parameters by 
optimization techniques was necessary. 

In future studies, by introducing an informational delay 
(Draye et al., 1997) or an artificial distance based on a Gaussian 
factor modulating the weights between the different neurons 
(Draye et al., 2002), it will be possible to analyze the self-tailored 
organization of the links between neurons and the possible emer- 
gence of specific topologies. In this case it will also be possible to 
select different modular architectures of the DRNN. 
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